#include <matrix.hpp>

#include <iostream>
#include <algorithm>

#include <iomanip>


#include <boost/lambda/lambda.hpp>

#include <random.hpp>



using namespace std;
using namespace sm;

int main()
{
	cout << setprecision(10);

	feng::variate_generator<long double> Ran;

	static_matrix<long double, 10, 10> m1;
	//generate(m1.rbegin(), m1.rend(), Ran());
	for_each( m1.begin(), m1.end(), boost::lambda::_1=Ran() );
	fill(m1.diag_rbegin(), m1.diag_rend(), -10.0 );
	fill(m1.upper_diag_rbegin(5), m1.upper_diag_rend(5), -1.0 );
	fill(m1.col_rbegin(3), m1.col_rend(3), 9.8);
	fill(m1.row_rbegin(7), m1.row_rend(7), -2.9);
	cout << "\nm1 = \n" << m1 << endl;

	static_matrix<long double, 11, 9> m2;
	fill(m2.rbegin(), m2.rend(), 10.0);
	cout << "\nm2 = \n" << m2 << endl;
	
	static_matrix<long double, 10, 10> m3;
	fill(m3.rbegin(), m3.rend(), Ran());
	fill(m3.diag_rbegin(), m3.diag_rend(), -9.0);
	fill(m3.upper_diag_rbegin(5), m3.upper_diag_rend(5), 17.0 );
	fill(m3.lower_diag_rbegin(1), m3.lower_diag_rend(1), 117.0 );
	fill(m3.col_rbegin(2), m3.col_rend(2), 7.8);
	fill(m3.row_rbegin(9), m3.row_rend(9), -1.3);
	cout << "\nm3 = \n" << m3 << endl;
	cout << "\n1/m3 = \n" << 1.0/m3 << endl;
	
	static_matrix<long double, 10, 10> m4;
	m4 = m3*m1;
	cout << "\nm4 = \n" << m4 << endl;

	static_matrix<long double, 10, 10> m5;
	m5 = m4/m1;
	cout << "\nm5 = \n" << m5 << endl;
	cout << "\nm5.det = \n" << m5.det() << endl;
	cout << "\nm3 = \n" << m3 << endl;
	cout << "\nm3.det = \n" << m3.det() << endl;
	
	static_matrix<long double, 10, 10> m5m5 = m5 * m5.invert();
	static_matrix<long double, 10, 10> m3m3 = m3 * m3.invert();

	cout << "\nm5*m5^{-1} = \n" << m5m5 << endl;
	cout << "\nm3*m3^{-1} = \n" << m3m3 << endl;

	cout << "\nm5*m5^{-1}.det = \n" << m5m5.det() << endl;
	cout << "\nm3*m3^{-1}.det = \n" << m3m3.det() << endl;

	/*
	static_matrix<long double, 500, 500> m6;
	//generate( m6.begin(), m6.end(), Ran() );
	for_each( m6.begin(), m6.end(), boost::lambda::_1=ran() );
	
	static_matrix<long double, 500, 500> m7 = m6.invert();
	static_matrix<long double, 500, 500> m8 = m6 * m7;
	cout << "\nUnit m8 = \n" << m8 << endl;
*/




	return 0;
}
